#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#define N 1024
//define struct
typedef struct matrix
{
    int rows;
    int **body;
}matrix;
//apply for space
matrix *malloc_matrix(int rows)
{
    int i;
	matrix *result;
    result=(matrix *)malloc(sizeof(matrix));
    result->body=(int **)malloc(sizeof(int)*rows);
    for(i=0;i<rows;i++)
        result->body[i]=(int *)malloc(sizeof(int)* rows);    
    result->rows=rows;
    return(result);
}
//release space
void free_matrix(matrix *m)
{
    int i;
    for(i=0;i<m->rows;i++)
    {
        free(m->body[i]);
    }
    free(m->body);
    free(m);
}
//addition
matrix *matrix_add(matrix *A,matrix*B)
{
    int i,j;
	matrix *result;
    result=malloc_matrix(A->rows);
    for(i=0;i<A->rows;i++)
        for(j=0;j<A->rows;j++)
            result->body[i][j]=A->body[i][j]+B->body[i][j];
    return(result);
}
//subtract
matrix *matrix_sub(matrix *A,matrix*B)
{
    int i,j;
	matrix *result;
    result=malloc_matrix(A->rows);
    for(i=0;i<A->rows;i++)
    {
        for(j=0;j<A->rows;j++)
        {
            result->body[i][j]=A->body[i][j]-B->body[i][j];

        }
    }
    return(result);
}
//divide origional matrix into four
matrix ** matrix_split(matrix *A)
{
    int i,j;
	matrix **ms;
	int rows;
    ms=(matrix **)malloc(sizeof(matrix*)*4);
    rows=A->rows/2;
    ms[0]=malloc_matrix(rows);
    ms[1]=malloc_matrix(rows);
    ms[2]=malloc_matrix(rows);
    ms[3]=malloc_matrix(rows);
    for(i=0;i<rows;i++)
        for(j=0;j<rows;j++)
        {
            ms[0]->body[i][j]=A->body[i][j];
			ms[1]->body[i][j]=A->body[i][j+rows];
			ms[2]->body[i][j]=A->body[i+rows][j];
			ms[3]->body[i][j]=A->body[i+rows][j+rows];
        }
    return(ms);
}
//splice matrix
matrix *matrix_merge(matrix *ms1,matrix *ms2,matrix *ms3,matrix *ms4)
{
    int i,j;
	matrix *result;
	result=malloc_matrix(ms1->rows+ms1->rows);
    for(i=0;i<ms1->rows;i++)
        for(j=0;j<ms1->rows;j++)
        {
            result->body[i][j]=ms1->body[i][j];
			result->body[i][j+ms1->rows]=ms2->body[i][j];
			result->body[i+ms1->rows][j]=ms3->body[i][j];
			result->body[i+ms1->rows][j+ms1->rows]=ms4->body[i][j];
        }
    return(result); 
}
//mult of second order matrix 
matrix *stramul(matrix *A,matrix *B,matrix *c)
{
	int pi1,pi2,pi3,pi4,pi5,pi6,pi7;
	pi1=(A->body[0][1]-A->body[1][1])*(B->body[1][0]+B->body[1][1]);
	pi2=(A->body[0][0]+A->body[1][1])*(B->body[0][0]+B->body[1][1]);
	pi3=(A->body[1][0]-A->body[0][0])*(B->body[0][0]+B->body[0][1]);
	pi4=(A->body[0][0]+A->body[0][1])*B->body[1][1];
	pi5=A->body[0][0]*(B->body[0][1]-B->body[1][1]);
	pi6=(B->body[1][0]-B->body[0][0])*A->body[1][1];
	pi7=(A->body[1][0]+A->body[1][1])*B->body[0][0];
	c->body[0][0]=pi1+pi2+pi6-pi4;
	c->body[0][1]=pi4+pi5;
	c->body[1][0]=pi6+pi7;
	c->body[1][1]=pi2+pi3+pi5-pi7;
	return(c);
}
//partitioning of matrix
matrix *matrixmul(matrix *A,matrix *B,matrix *c)
{	
	matrix *matrix_mul(matrix *A,matrix *B);
	matrix *p1,*p2,*p3,*p4,*p5,*p6,*p7;
	matrix *c11,*c12,*c21,*c22;
	matrix **sa,**sb;
	int rows;
	rows=A->rows/2;
	c11=malloc_matrix(rows);
	c12=malloc_matrix(rows);
	c21=malloc_matrix(rows);
	c22=malloc_matrix(rows);
	sa=matrix_split(A);sb=matrix_split(B);
	p1=matrix_mul(sa[0],matrix_sub(sb[1],sb[3]));
	p2=matrix_mul(matrix_add(sa[0],sa[1]),sb[3]);
	p3=matrix_mul(matrix_add(sa[2],sa[3]),sb[0]);
	p4=matrix_mul(sa[3],matrix_sub(sb[2],sb[0]));
	p5=matrix_mul(matrix_add(sa[0],sa[3]),matrix_add(sb[0],sb[3]));
	p6=matrix_mul(matrix_sub(sa[1],sa[3]),matrix_add(sb[2],sb[3]));
	p7=matrix_mul(matrix_sub(sa[0],sa[2]),matrix_add(sb[0],sb[1]));
	c11=matrix_add(matrix_add(p5,p4),matrix_sub(p6,p2));
	c12=matrix_add(p1,p2);
	c21=matrix_add(p3,p4);
	c22=matrix_add(matrix_sub(p5,p3),matrix_sub(p1,p7));
	c=matrix_merge(c11,c12,c21,c22);
	free(p1);
	free(p2);
	free(p3);
	free(p4);
	free(p5);
	free(p6);
	free(p7);
	free(c11);
	free(c12);
	free(c21);
	free(c22);	
    return(c);
}
//traditional matrix mult
matrix *matrix_mul(matrix *A,matrix *B)
{
	matrix *c;
	c=malloc_matrix(A->rows);
	if(A->rows==2)
		c=stramul(A,B,c);
	else
		c=matrixmul(A,B,c);
	return(c);
}
//print matrix 
void print_matrix(matrix *A)
{
    int i,j;
    for( i=0;i<A->rows;i++)
    {
        for( j=0;j<A->rows;j++)
        {
            printf("%d\t",A->body[i][j]);	

        }
        printf("\n");
    }
}
matrix *ordmult(matrix *A,matrix *B)
{
	int i,j,k;
	matrix *c;
	c=malloc_matrix(A->rows);
	for(i=0;i<A->rows;i++)
		for(j=0;j<B->rows;j++){
			c->body[i][j]=0;
				for(k=0;k<A->rows;k++)
					c->body[i][j]=c->body[i][j]+A->body[i][k]*B->body[k][j];
		}
	return(c);
}
void main()
{
    matrix *A,*B,*c1,*c2;
    int i,j;
	long double dura1,dura2;
	clock_t start,finish;
	A=malloc_matrix(N);B=malloc_matrix(N);
	srand((unsigned int)time(NULL));
    for( i=0;i<A->rows;i++)
    {
        for(j=0;j<A->rows;j++)
        {
            A->body[i][j]=rand()%10;
        }
    } 
    for( i=0;i<B->rows;i++)
    {
        for(j=0;j<B->rows;j++)
        {
            B->body[i][j]=rand()%10;
        }
    }
	start=clock();
    c1=matrix_mul(A,B);
	finish=clock();
	dura1=(double)(finish-start)/CLOCKS_PER_SEC;
	start=clock();
	c2=ordmult(A,B);
	finish=clock();
	dura2=(double)(finish-start)/CLOCKS_PER_SEC;
//	print_matrix(c2);
//	print_matrix(c1);
	printf("%lf\n",dura1);
	printf("%lf\n",dura2);
	system("pause");
}
